arXiv:1505.06968vl [cond-mat.soft] 24 May 2015 


1 


A Nonlinear Boundary Condition for Continuum Models of 

Biomolecular Electrostatics 

J. P. Bardhan^ (corresponding and presenting author), D. A. Tejani^, N. S. Wieckowski^, A. 

Ramaswamy^, and M. G. Knepley^ 

Northeastern University, Boston, MA, USA (j .bardhcinOneu.edu) 

^University of Chicago, Chicago, IL, USA (knepley@ci.uchicago.edu) 


Abstract — Understanding the behavior of biomolecules such as proteins requires understand¬ 
ing the critical influence of the surrounding fluid (solvent) environment—water with mobile salt 
ions such as sodium. Unfortunately, for many studies, fully atomistic simulations of biomolecules, 
surrounded by thousands of water molecules and ions are too computationally slow. Continuum- 
solvent models based on macroscopic dielectric theory (e.g. the Poisson equation) are popular 
alternatives, but their simplicity fails to capture well-known phenomena of functional signifi¬ 
cance. For example, standard theories predict that electrostatic response is symmetric with 
respect to the sign of an atomic charge, even though response is in fact strongly asymmetric if 
the charge is near the biomolecule surface. In this work, we present an asymmetric continuum 
theory that captures the essential physical mechanism—the finite size of solvent atoms—using a 
nonlinear boundary condition (NLBC) at the dielectric interface between the biomolecule and sol¬ 
vent. Numerical calculations using boundary-integral methods demonstrate that the new NLBC 
model reproduces a wide range of results computed by more realistic, and expensive, all-atom 
molecular-dynamics (MD) simulations in explicit water. We discuss model extensions such as 
modeling dilute-electrolyte solvents with Debye-Huckel theory (the linearized Poisson-Boltzmann 
equation) and opportunities for the electromagnetics community to contribute to research in this 
important area of molecular nanoscience and engineering. 


1. INTRODUCTION 

Protein structure and function are determined in part by electrostatic interactions between the 
protein’s atomic charges and the surrounding biological fluid (solvent), a complex mixture of polar 
water molecules and dissolved charges such as sodium and potassium ions [T]. For over a century, 
biological scientists have modeled these interactions using macroscopic continuum models based on 
the Poisson-Boltzmann equation or Debye-Hiickel theory mm- These popular mean-field theories 
assume that solvent molecules are infinitely small compared to the biomolecule solute mil], a drastic 
simplification critical to enable pioneering theoretical studies using spherical and ellipsoidal models 
of protein shapes [5]- However, its justifications are increasingly questionable for atomistically 
detailed protein structures and predictions of binding affinities (binding free energies), which are 
of enormous value to understanding the molecular basis of disease and for designing improved 
therapeutic drugs. 

In contrast, exponential growth in computing capabilities has enabled large-scale molecular- 
dynamics (MD) simulations that model the surrounding solvent in fully atomistic, explicit de¬ 
tail [6]. These calculations can provide far more realistic insights, but unfortunately require hun¬ 
dreds or thousands of times the computational resources that continuum models need. A rigorous 
statistical-mechanical argument establishes the connection between atomistic and continuum mod¬ 
els mu], creating an opportunity to develop more accurate continuum models by comparison to 
MD simulations, where new continuum theories can be tested against a wide array of challenging 
“computational experiments” which are unrealizable in real-world laboratories [niHiiiniiii]. 

We have been advancing multiple approaches to improve continuum models [Hiiiiiiiaiii]. One 
approach replaces the traditional model of the solvent as a macroscopic dielectric continuum (the 
familiar relation D(r) = e{r)E{r)), and instead models it as a nonlocal dielectric material [151 HHTB]. 
which limits short-range dielectric response using a convolution D{r) = “ r'\)E{r')dV'. 

On the other hand, nonlocal models, and even many sophisticated nonlinear models, are still 
symmetric with respect to the sign of the protein charges. In reality, however, negative and positive 
charges near the solute surface produce substantially different responses. Asymmetry results from 
both water structure at the interface, and the fact that water hydrogens are much smaller than 
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water oxygens [m El 11 ng HI [201 HI]. We have demonstrated, using atomistic MD calculations, 
that charge-sign asymmetries can be accurately reproduced if the electric potential is modeled as a 
piecewise-linear plus constant (i.e. piecewise-affine) function of the charge |11] . The constant term 
in this model represents the interface potential induced by water structure around a completely 
uncharged version of the solute [SUES], and the discontinuity in response coefficient occurs where 
the solute charge changes sign m- 

The few theories that directly address hydration asymmetry [231 [T^ EQ] are not actually Poisson 
models, but generalizations of the Born-ion problem (the spherically symmetric case of a sphere 
with a central charge). Recently we proposed the first successful asymmetric Poisson theory that 
can be solved for complex molecular geometries m, by translating the existing models’ physical 
insights into a boundary-integral equation (BIE) formulation of the Poisson problem [241125j . This 
led to a modified BIE formulation in which we replaced the usual Maxwell boundary condition for 
the continuity of the normal flux with a nonlinear boundary condition (NLBC) [14] . Calculations 
showed that the new model successfully reproduces highly accurate MD calculations on a wide 
range of challenging test cases m- 

In the present paper, we improve the new NLBC model by ensuring that the potential outside 
the solute protein satisfies Gauss’s law sufficiently far from the solute. To demonstrate the new 
model, we have implemented a boundary-element method solver in MATLAB (all of the source code 
and data required to reproduce the figures in this paper are freely and publicly available online [26]). 
Calculations demonstrate that the improved model is more accurate for monovalent atomic ions, and 
illustrate that charge-sign hydration asymmetry effects are substantial for surface charges, with the 
magnitude of asymmetry decreasing rapidly for charges further from the solute-solvent interface. In 
particular, we predict that for surface charges in a sphere, the energetic difference between positive 
and negative charges does not depend strongly on the size of the sphere; this result implies that 
asymmetry is an essential piece of physics to include in models of molecular electrostatics. The 
next section describes the traditional, charge-sign symmetric continuum model, our modifications to 
incorporate asymmetric response, and a boundary-integral approach for solving the NLBC model. 
Section El describes a MATLAB implementation and presents computational results for single-atom 
ions as well as large spheres approximating full-sized proteins. The paper concludes in Section Ej 
with a summary and discussion. 

2. THEORY 

2.1. Symmetric Continuum Model and Boundary-Integral Method 

We consider a protein in an infinite dilute electrolyte solution. In the exterior region, which we label 
/, the potential obeys the linearized Poisson-Boltzmann equation where k represents 

the inverse Debye screening length [T], and the dielectric constant is labeled e/ (often taken to be 80, 
approximately that of bulk water). A thin shell of ion-free solvent separates the electrolyte solution 
from the protein. This region II, labeled the Stern or ion-exclusion layer, is a few Angstroms 
in width and simplistically models the finite size of the ions dissolved in the electrolyte. Here 
the potential obeys the Laplace equation and the dielectric constant is e// (usually e// = e/). The 
protein interior, labeled III, is a low-dielectric medium (e/// is usually between 2 and 8) containing 
N discrete point charges, and the potential satisfies a Poisson equation V‘^(p{r) = — 
where qi and r* specify the ith charge. The boundary a separates the protein region III and the 
Stern layer II, and the boundary b separates the Stern layer from the electrolyte solvent I. The 
potential is assumed to decay to zero suitably fast as r ^ oo, and the boundary conditions are 


P’lllira) = P>Il{ra) 

( 1 ) 

d^piiiira) dip 11 {r a) 

r. = eil 

on on 

( 2 ) 


( 3 ) 

difni^b) dipi{rb) 

dn dn ■ 

( 4 ) 


The flux conditions Eqs. [2] and 0] are what we call the standard Maxwell boundary conditions 
(SMBC). By applying Green’s theorem in these three regions and taking suitable limits as the field 
point approaches each region’s bounding surface or surfaces, we obtain the coupled BIE system m 
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[Ml [28]: 
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where Gx,ij and Kx,ij represent the single- and double-layer boundary-integral operators associated 
with the Green’s function of region X that map from a distribution on boundary j to the potential 
on boundary i. 

2.2. Renormalized Nonlinear Boundary Condition Model 

Our original work on the NLBC model employed only a single interface, the protein solvent-solute 
boundary a. For consistency, the regions it separates are still labeled III (solute) and II (solvent), 
and instead of Eq. [2] we have the nonlinear flux boundary condition 




where the field-dependent nonlinear function / is 
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The first term of Eq. [7] represents the SMBC, and the last term ensures that the NLBC model 
recovers the standard model for weak electric fields, i.e. as —>■ 0 . As shown in our earlier work, 
the NLBC has only three free parameters a, /3, and 7 , which can be parameterized robustly against 
atomistic calculations. Numerical simulations using the NLBC model showed excellent agreement 
with atomistic results [Il|. However, outside the solute, the potentials generated by this model 
fail to satisfy Gauss’s law, as can be seen by considering a Born ion, i.e. a sphere with central 
charge. In particular, solutions satisfying SMBC automatically satisfy Gauss’s law, whereas the 
NLBC cannot simultaneously satisfy Gauss’s law and provide asymmetric response. This dehciency 
of central importance for problems involving multiple solute molecules, e.g. protein-drug binding; 
we propose to solve it by including a compensating charge distribution a few Angstroms away from 
the solute-solvent interface, at the Stern layer. This compensating charge ensures that Gauss’s law 
is satished outside the second boundary. The modihed boundary condition at b ensures that the 
model recovers the expected macroscopic dielectric behavior far from the protein surface, and is 
written 
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Defining si and S 2 so that Eq. [ 8 ] can be written si ~ '^ 2 ^, the system of coupled BIEs is 
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Using Green’s theorem again to determine the reaction potential yj^“‘’(g) in the protein due to the 
solvent, the electrostatic charging free energy is then calculated as AG = E) 

where we have modeled the static (interface) potential as a constant, see e.g. [S] [Ml [II] (following 
our previous work, we model = 10.7 kcal/mol/e). In contrast, the standard Poisson theory 

gives an energy AG = Lq where L is a symmetric negative dehnite operator. 


3. COMPUTATIONAL RESULTS 
3.1. Numerical Implementation 

The full MATLAB code to reproduce the calculations and figures in this paper are freely available 
online at http://www.bitbucket.org/jbardhan/piersl5-code, Our boundary-element method 
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from earlier work was extended from using only point-based discretizations of the relevant bound¬ 
aries and unknown surface distributions, to use triangular boundary elements with piecewise- 
constant basis functions and centroid collocation. Following our earlier work, we use Picard it¬ 
eration to solve the nonlinear BIE problem. Calculations on spherical molecules used the earlier 
point-based implementation for verification against earlier results. Calculations for ellipsoids used 
triangular meshes derived from the mesh obtained from MATLAB’s ellipsoid command, which 
takes as input the ellipsoid semi-axis lengths and a number n of subdivisions, and returns three 
(n-|- 1 ) by (n-|- 1 ) matrices (for the x, y, and 2 : coordinates of the mesh vertices), which define planar 
quadrilaterals and triangles by subdividing the ellipsoid in angular coordinates (lines of latitude 
and longitude). By iterating over the quadrilaterals and subdividing, we obtain a triangular mesh 
suitable for our existing MATLAB implementation of the boundary-element method. 

The ellipsoidal shape approximations are calculated using standard methods [29l |30] : the 
molecule, e.g. a protein, is dehned as a set of atoms which are dehned as spheres at specihed 
locations and with specified radii, and with each possessing a single point charge at its center. The 
N atomic positions are represented by rj = [xi‘,yi]Zi] and the radii are [oi, 02 ,...., oat]. In order 
to estimate protein shape as an ellipsoid, we translate the molecule so its center of mass is at the 
origin, where the “mass” of the molecule is modeled as proportional to the sum of the atomic 
volumes, i.e. M = center of mass is dehned as Vc = (Af“^) ^^;^(a?)(rj), and so 

we translate the atoms according to r( = — Tc. Dropping the prime and referring only to the 

translated coordinates, the components of the molecule’s inertia tensor I are calculated as 


hi = '^{mi{yf + zf + ^a^)) ( 10 ) 

i=l 

h 2 = ^{rriiixj + zf + ^aj)) ( 11 ) 

i=l 

hs = ^{rriiix'f + yf + ^a^)) (12) 

i=l 

lu = '^{mxiyi) (13) 

i=l 

Ii3 = '^{rriiXiZi) (14) 

i=l 

h3 = '^{rriiViZi), (15) 

i=l 


and by symmetry /21 = / 12 , ^31 = / 13 , and 1^2 = hs- The principal moments of inertia of the 
molecule are the eigenvalues of I. Choosing Ixx, lyy, and Izz such that Ixx < lyy < hz, we 
hnd the semi-axes a, b, and c of an ellipsoid that has the same weight M and principal moments 

of inertia [29l |30j. This leads to a = ^ = \/ ^(hx - lyy + hz), and 

c = This ellipsoid is the simplest anisotropic approximation to the shape of 

the bio-molecule under consideration. 

The results in Figure [T](a) represent continuum-model calculations of monovalent ions (-|-le and 
— le charges) as a function of ion radius. For comparison, reference data for biologically relevant 
ions obtained from all-atom molecular dynamics (MD) calculations [T3] are plotted as symbols. The 
renormalized NLBC (thick curves) clearly fit the MD results better than the original NLBC model 
(thin curves), even though no new parameters have been introduced. The symmetric Poisson model 
is substantially incorrect; the most frequent approaches for ions involve adjusting atomic radii on 
an atom by atom basis, but this is not possible for multi-atom solutes [T9l[20]. Figure dKb) plots 
the deviations between the new NLBC model and a Poisson model that involves sign-symmetric 
dielectric response, but incorporates the static potential contribution: that is, the results plot the 
effects of asymmetry in the dielectric response at the molecule surface. The sample problems in this 
figure are a surface charge in a sphere of varying radius (1.5 A from the surface), or a buried charge 
at its center. As expected, the buried charge experiences essentially symmetric response once it 
is more than a few Angstroms from the interface. In contrast, the surface charge experiences a 
surprisingly large asymmetry that is essentially constant even as the sphere radius increases. The 
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magnitude of this asymmetry is in line with previous MD calculations and the persistence of 
this large deviation from standard models, even for protein-sized molecules, suggests that including 
accurate models of asymmetry should be a main goal in developing improved continuum theory. 
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Figure 1: Charging free energies for single charges in spheres, (a) Common monovalent ions, calculated 
using all-atom molecular dynamics (symbols) [M] , standard Poisson theory (black dashed line), the original 
NLBC model [14] (thin lines), and the new NLBC model presented in this work (thick lines), (b) Charging 
free energy for a single charge in a sphere, as a function of sphere radius, charge value, and charge location. 
Buried charges are situated at the sphere center; surface charges are located 1.5 A from the sphere surface. 


4. DISCUSSION 


In this work we have presented a continuum-electrostatic model for molecular solvation that models 
biologically important hydration asymmetry-a fundamentally non-continuum phenomenon-using 
an effective nonlinear boundary condition (NLBC). The NLBC replaces the traditional Maxwell 
boundary condition that enforces continuity of the electric flux at the interface, and the present 
model improves on our original work m by ensuring that Gauss’s law holds outside of the first 
shell of solvent molecules surrounding molecular solutes. The use of an effective boundary condi¬ 
tion at the molecule-solvent interface represents a new frontier in biomolecular modeling, and was 
inspired by a long history of approximate boundary conditions in electromagnetic theory, e.g. m 
and rarefied gases [SS] |33] . We hope that the present work will encourage experts in the electro¬ 
magnetics community to contribute their insights and experience to support deeper understanding 
of molecular electrostatics, whether in improving models, analyzing their implications, or solving 
them numerically. 

Our NLBC model is the first asymmetric theory that uses actual Poisson continuum theory, and 
in the simple case of a single ion, reproduces empirical and semi-empirical theories developed over 
decades of research into electrostatic asymmetry nziiisiEn]. We note that the present model can 
treat dilute electrolyte solutions using the linearized Poisson-Boltzmann equation BE] simply by 


modeling the Green’s function of the outermost region with the LPBE Green’s function ^ 

instead of the of the Laplace Green’s function. To demonstrate the new model’s accuracy, 

we have calculated the electrostatic charging free energies for single-atom ions and compared our 
results to more accurate, and much more computationally demanding, atomistic molecular dynam¬ 
ics calculations that include thousands of explicit water molecules. The results for ions in Figure 1 
indicate that the new model exhibits substantially better accuracy than the original, and the en¬ 
forcement of Gauss’s law is in fact even more important when modeling electrolyte solutions using 
the LPBE (results not shown). Galculations of amino-acid charging free energies, and the differ¬ 
ences due to protonation or deprotonation, illustrate that the magnitude of asymmetric response 
decays rapidly with an atomic charge’s distance from the solute-solvent interface UDEli. Our 
numerical calculations employed boundary-integral methods with simple model geometries, such as 
spheres for the monatomic ions and ellipsoids to model the amino acids. These ellipsoids represent 
simple shape approximations [29l |30] and we expect that they will be useful for fast approximate 
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calculations such as in implicit-solvent molecular dynamics [Mll25l|35l[36]. Calculations for atom¬ 
istic models of large molecules such as proteins will require fast, parallel boundary-element method 
solvers [24t 128] , and implementation of such software represents an area of ongoing work. 
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